Numerical integration of orbital motion

Example1: Orbit around Earth (central force)

Here is a simple example of how to use the force functions (plus a few other ones) defined in Celestlab to integrate the motion.

Only the central force is considered. The trajectory is compared with a purely Keplerian trajectory.

// -----------------------------------------------------
// Utilities
// -----------------------------------------------------
// State evolution: dy/dt = f(t, y). 
// y = [pos; vel] (state vector)
// t in seconds
// model = central force
function [ydot]=fct(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 
  ydot = zeros(y);
  ydot(1:3,:) = vel;
  ydot(4:6,:) = CL_fo_centralAcc(pos); 
endfunction

// -----------------------------------------------------
// Initializations, integration
// -----------------------------------------------------
// initial state
t0 = 0;  
pos0 = [7000.e3; 0; 0]; 
vel0 = [0; 5.e3; 5.e3];
y0 = [pos0; vel0]; 

t = (0:180:10000); // time instants (seconds)

// integration
rtol = 1.e-12 * [1;1;1;1;1;1];
atol = 1.e-6 * [1;1;1;1.e-3;1.e-3;1.e-3]; 

y = ode(y0, t0, t, rtol, atol, fct);
pos = y(1:3,:); 
vel = y(4:6,:); 

// Comparison with Keplerian motion 
kep0 = CL_oe_car2kep(pos0, vel0); 
kep = CL_ex_kepler(0, kep0, t/86400); 
[pos_k, vel_k] = CL_oe_kep2car(kep); 

// -----------------------------------------------------
// Plots
// -----------------------------------------------------

scf();
plot(t/86400, CL_norm(pos-pos_k),"k"); 
xtitle("Error on position (m)", "Time (days)"); 
CL_g_stdaxes()

scf();
plot(t/86400, CL_norm(vel-vel_k),"k"); 
xtitle("Error on velocity (m/s)", "Time (days)"); 
CL_g_stdaxes()

Example2: Orbit around Earth (central force + J2..J6)

Here is another example of integration of orbital motion.

Initial position and velocity are initialized from mean orbital elements (Eckstein-Hechler model). The trajectory is compared with that obtained by Eckstein-Hechler. Energy is also computed, which gives an idea of the integration accuracy.

// -----------------------------------------------------
// Utilities
// -----------------------------------------------------
// State evolution: dy/dt = f(t, y). 
// y = [pos; vel] (state vector)
// t in seconds
// model = central force + zonal terms (J2:J6)
function [ydot]=fct(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 
  ydot = zeros(y);
  ydot(1:3,:) = vel;
  ydot(4:6,:) = CL_fo_centralAcc(pos) + CL_fo_zonHarmAcc(pos, 2:6); 
endfunction

// Potential (for integration checking)
function [pot]=potential(pos) 
  pot = CL_fo_centralPot(pos) + CL_fo_zonHarmPot(pos, 2:6); 
endfunction

// -----------------------------------------------------
// Initializations, integration
// -----------------------------------------------------
// orbital elements = sma, ex, ey, inc, raan, psom 
t0 = 0; 
param0_mean = [7000.e3; 0; 0.001; 98*%pi/180; 0; 0]; 
param0_osc = CL_ex_propagate("eckhech", "cir", t0, param0_mean, t0, "o"); 
[pos0, vel0] = CL_oe_cir2car(param0_osc); 
y0 = [pos0; vel0]; // state vector
t = (0:180:86400); // time instants (seconds)

// integration
rtol = 1.e-12 * [1;1;1;1;1;1];
atol = 1.e-6 * [1;1;1;1.e-3;1.e-3;1.e-3]; 

y = ode(y0, t0, t, rtol, atol, fct);
pos = y(1:3,:); 
vel = y(4:6,:); 

// Conversion to mean orbital elements 
param_osc = CL_oe_car2cir(pos, vel); 
param_mean = CL_ex_osc2mean("eckhech", "cir", param_osc); 

// propagation with Eckstein-Hechler 
[param_mean_eh, param_osc_eh] = CL_ex_propagate("eckhech", "cir", 0, ..
                                param0_mean, t/86400, "mo"); 
[pos_eh, vel_eh] = CL_oe_cir2car(param_osc_eh); 

// Energy
mu = CL_dataGet("mu"); 
energy = CL_dot(vel)/2 - potential(pos);

// -----------------------------------------------------
// Plots
// -----------------------------------------------------

scf();
plot(t/86400, CL_norm(pos-pos_eh),"k"); 
xtitle("Position difference with Eckstein-Hechler J2..J6 (m)", "Time (days)"); 
CL_g_stdaxes()

scf();
plot(param_mean(2,:), param_mean(3,:), "b"); 
plot(param_mean_eh(2,:), param_mean_eh(3,:), "or"); 
xtitle("Mean eccentricity vector"); 
CL_g_stdaxes()

scf();
plot(t/86400, -mu./(2*energy) + mu./(2*energy(1)),"k"); 
xtitle("Energy variations (m)", "Time (days)"); 
CL_g_stdaxes()

Example3: Orbit around Earth (all harmonics)

Here is another example of integration of orbital motion.

The trajectory is integrated in the inertial frame and in the rotating frame. The two solutions are compared with each other. Energy is also computed, which gives an idea of the integration accuracy.

// -----------------------------------------------------
// Utilities
// -----------------------------------------------------
// State evolution: dy/dt = f(t, y). 
// y = [pos; vel] (state vector)
// t in seconds
// model = central force + spherical harmonics (10x10)

// Case 1: 
// y = [pos; vel] relative to inertial frame (components in inertial frame)
// Inertial frame = ECF frozen at T0
// => Conversion from ECI to ECF required
function [ydot]=fct1(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 
  ydot = zeros(y);
  ydot(1:3,:) = vel;
  mat = CL_rot_angles2matrix(3, OMEGA(3)*t); 
  pos_ecf = mat * pos; 
  acc_sphharm = mat' * CL_fo_sphHarmAcc(pos_ecf, [10,10]);
  ydot(4:6,:) = CL_fo_centralAcc(pos) + acc_sphharm; 
endfunction

// Case 2: 
// y = [pos; vel] relative to rotating frame (components in rotating frame)
// => Apparent acceleration added
function [ydot]=fct2(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 
  ydot = zeros(y);
  ydot(1:3,:) = vel;
  acc_sphharm = CL_fo_sphHarmAcc(pos, [10,10]);
  ydot(4:6,:) = CL_fo_centralAcc(pos) + acc_sphharm + ..
                CL_fo_apparentAcc(pos, vel, OMEGA); 
endfunction

// Potential (for integration checking)
// pos: relative to rotating frame
function [pot]=potential(pos) 
  pot = CL_fo_centralPot(pos) + CL_fo_sphHarmPot(pos, [10,10]) + ..
        CL_fo_centrifugPot(pos, OMEGA); 
endfunction

// -----------------------------------------------------
// Initializations, integration
// -----------------------------------------------------

// Earth angular velocity (ECI -> ECF)
// NB: value could be also be obtained by CL_mod_siderealTime()
// Components can be considered as given in ECF or ECI
OMEGA = [0; 0; CL_dataGet("rotrBody")];  

// initial date (days from 1950.0) 
T0 = 20000;  

// initial state
t0 = 0;  
pos0_eci = [42166.e3; 0; 0]; 
vel0_eci = [0; 3075; 0];

t = (0:300:2*86400); // time instants (seconds)

// integration
rtol = 1.e-12 * [1;1;1;1;1;1];
atol = 1.e-6 * [1;1;1;1.e-3;1.e-3;1.e-3]; 

// ------ CASE 1 ------
// Integration in inertial frame (=ECF frozen at T0)
disp("Integration - case1"); 
mat = CL_fr_convertMat("ECI", "ECF", T0); 
y0 = [mat * pos0_eci; mat * vel0_eci]; 

y = ode(y0, t0, t, rtol, atol, fct1);

mat = CL_rot_angles2matrix(3, OMEGA(3)*t); 
[pos_ecf1, vel_ecf1] = CL_rot_pvConvert(y(1:3,:), y(4:6,:), mat, OMEGA*ones(t)); 

// potential   
[pot1] = potential(pos_ecf1); 
energy1 = CL_dot(vel_ecf1)/2 - pot1; 

// ------ CASE 2 ------
// Integration in rotating frame
disp("Integration - case2"); 
[pos0, vel0] = CL_fr_convert("ECI", "ECF", T0, pos0_eci, vel0_eci); 
y0 = [pos0; vel0]; 

y = ode(y0, t0, t, rtol, atol, fct2);
pos_ecf2 = y(1:3,:); 
vel_ecf2 = y(4:6,:); 

// potential   
[pot2]=potential(pos_ecf2); 
energy2 = CL_dot(vel_ecf2)/2 - pot2; 

// -----------------------------------------------------
// Plots
// -----------------------------------------------------
mu = CL_dataGet("mu"); 

scf();
plot(t/86400, energy1 - energy1(1),"b"); 
plot(t/86400, energy2 - energy2(1),"r"); 
xtitle("Energy variations (m2/s2)", "Time (days)"); 
CL_g_stdaxes(); 

scf();
plot(t/86400, CL_norm(pos_ecf2-pos_ecf1), "k"); 
xtitle("Position difference (m)", "Time (days)"); 
CL_g_stdaxes();

Example4: Orbit around Lagrange point (Halo orbit)

Here is another example of integration of orbital motion.

A Halo orbit in computed using the dedicated function, and is compared with the solution of a numerical integration (using the standard force functions). We observe that the integrated trajectory diverges after about 1.5 orbit due to the great sensitivity around the Lagrange point.

// -----------------------------------------------------
// Utilities
// -----------------------------------------------------
// State evolution: dy/dt = f(t, y): 
// Acceleration in rotating frame (center = barycenter of body1/body2)
// NB: uses "global" variables (POS1, POS2...) 
function [ydot]=fct(t, y)
  pos = y(1:3,:); 
  vel = y(4:6,:); 

  ydot = zeros(y);
  ydot(1:3,:) = vel;

  n = size(y, 2); 
  acc = CL_fo_centralAcc(pos - repmat(POS1,1,n), MU1) + ..
        CL_fo_centralAcc(pos - repmat(POS2,1,n), MU2) + ..
        CL_fo_apparentAcc(pos, vel, [0;0;OMEGA]); 

  ydot(4:6,:) = acc; 
endfunction

// Energy: kinetic energy minus potentil in rotating frame
function [e]=energy(pos, vel)
  n = size(pos, 2); 

  pot = CL_fo_centralPot(pos - repmat(POS1,1,n), MU1) + ..
        CL_fo_centralPot(pos - repmat(POS2,1,n), MU2) + ..
        CL_fo_centrifugPot(pos, [0;0;OMEGA]); 

  e = CL_dot(vel)/2 - pot; 
endfunction

// Initialisation of Lagrange point
env = CL_3b_environment("S-EM", "L2");

// Global variables for integration
// NB: Normalized (adimensional) coordinates
ratio = env.MU; 
POS1 =  -ratio * [1;0;0] / (1 + ratio); // Position of body1
POS2 = [1;0;0] / (1 + ratio); // Position of body2
OMEGA = 1; 
MU1 = 1 / (1 + ratio); // MU of body1
MU2 = 1 * ratio / (1 + ratio); // MU of body2

// Computation of a Halo orbit
Az = 150.e6 / env.D; // adimensional
direction = 0; // North
t_orb = linspace(0,1.5*360*86400,200) * env.OMEGA; // "times" (adimensional)
[orb,omega] = CL_3b_halo(env,Az,direction,t_orb); // "pv" (adimensional)

// Numerical integration starting from initial "time"
t0 = t_orb(1); 
y0 = orb(:,1); 

rtol = 1.e-12 * ones(6,1); 
atol = 1.e-12 * ones(6,1); 

y = ode(y0, t0, t_orb, rtol, atol, fct);
pos = y(1:3,:); 
vel = y(4:6,:); 

// Plots of position - normalized values
scf(); 
for k = 1:3
 subplot(3,1,k); 
 plot(t_orb, orb(k,:), "r"); 
 plot(t_orb, pos(k,:), "b"); 
 xtitle("Coord " + string(k)); 
end

// Plot of energy - normalized values
e = energy(pos, vel); 
scf(); 
xtitle("Energy"); 
plot(t_orb, e); 

// Plot of 3D trajectory - normalized values
scf(); 
xtitle("3D trajectory"); 
param3d(orb(1,:), orb(2,:), orb(3,:)); 
h = CL_g_select(gce()); 
h.foreground = 5; 
h.thickness = 2; 
param3d(pos(1,:), pos(2,:), pos(3,:)); 
h = CL_g_select(gce()); 
h.foreground = 2;

Report an issue
<< Frame conversions	Cookbook	Rotations >>